#ifndef _CUDA_TIJK_H_
#define _CUDA_TIJK_H_
#include <math.h>

#define _CUDA_TIJK_RANK1_eps_start	0.000000001f
#define _CUDA_TIJK_RANK1_eps_impr	0.0001f
#define _CUDA_TIJK_RANK1_beta		0.3f
#define _CUDA_TIJK_RANK1_gamma		0.9f
#define _CUDA_TIJK_RANK1_sigma		0.5f
#define _CUDA_TIJK_RANK1_maxtry		50

#define _CUDA_TIJK_RANKK_eps_res	0.000000001f
#define _CUDA_TIJK_RANKK_eps_impr	0.001f
#define _CUDA_TIJK_RANKK_pos		1

typedef struct cuda_tijk_type_t {
	unsigned int order; /* number of tensor indices */
	unsigned int dim; /* dimension of each axis (only square tensors supported) */
	unsigned int num; /* unique number of components */
#define CUDA_TIJK_TYPE_MAX_NUM 28
} cuda_tijk_type;

#define _CUDA_TIJK_6O3D_SYM_TSP(A, B)                                    \
	((A)[0]*(B)[0]+(A)[21]*(B)[21]+(A)[27]*(B)[27]+                   \
	6*((A)[1]*(B)[1]+(A)[2]*(B)[2]+(A)[15]*(B)[15]+(A)[20]*(B)[20]+  \
	(A)[22]*(B)[22]+(A)[26]*(B)[26])+                             \
	15*((A)[3]*(B)[3]+(A)[5]*(B)[5]+(A)[10]*(B)[10]+(A)[14]*(B)[14]+ \
	(A)[23]*(B)[23]+(A)[25]*(B)[25])+                            \
	30*((A)[4]*(B)[4]+(A)[16]*(B)[16]+(A)[19]*(B)[19])+              \
	20*((A)[6]*(B)[6]+(A)[9]*(B)[9]+(A)[24]*(B)[24])+                \
	60*((A)[7]*(B)[7]+(A)[8]*(B)[8]+(A)[11]*(B)[11]+                 \
	(A)[13]*(B)[13]+(A)[17]*(B)[17]+(A)[18]*(B)[18])+            \
	90*(A)[12]*(B)[12])                                              \

#define CUDA_TEN_T2M(m, t) ( \
   (m)[0] = (t)[1], (m)[1] = (t)[2], (m)[2] = (t)[3], \
   (m)[3] = (t)[2], (m)[4] = (t)[4], (m)[5] = (t)[5], \
   (m)[6] = (t)[3], (m)[7] = (t)[5], (m)[8] = (t)[6] )

/* ell */
#define CUDA_ELL_3V_SCALE(v2, a, v1) ((v2)[0] = (a)*(v1)[0], (v2)[1] = (a)*(v1)[1], (v2)[2] = (a)*(v1)[2])
#define CUDA_ELL_3V_DOT(v1, v2) ((v1)[0]*(v2)[0] + (v1)[1]*(v2)[1] + (v1)[2]*(v2)[2])
#define CUDA_ELL_3V_LEN(v) (sqrtf(CUDA_ELL_3V_DOT((v),(v))))
#define CUDA_ELL_3V_SCALE_ADD2(v2, s0, v0, s1, v1) \
  ((v2)[0] = (s0)*(v0)[0] + (s1)*(v1)[0],     \
   (v2)[1] = (s0)*(v0)[1] + (s1)*(v1)[1],     \
   (v2)[2] = (s0)*(v0)[2] + (s1)*(v1)[2])
#define CUDA_ELL_3V_NORM(v2, v1, length) (length = CUDA_ELL_3V_LEN(v1), CUDA_ELL_3V_SCALE(v2, 1.0f/length, v1))
#define CUDA_ELL_3V_COPY(v2, v1) ((v2)[0] = (v1)[0], (v2)[1] = (v1)[1], (v2)[2] = (v1)[2])
#define CUDA_ELL_3V_INCR(v2, v1) ((v2)[0] += (v1)[0], (v2)[1] += (v1)[1], (v2)[2] += (v1)[2])

#endif
